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Abstract 

New Physics that becomes relevant at some high scale A beyond the experimental reach, can 
be described in the effective theory approach by adding higher-dimensional operators to the 
Standard Model (SM) Lagrangian. In Higgs pair production through gluon fusion, which 
gives access to the trilinear Higgs self-coupling, this leads not only to modifications of the SM 
couplings but also induces novel couplings not present in the SM. For a proper prediction of 
the cross section, higher order QCD corrections that are important for this process, have to 
be taken into account. The various higher-dimensional contributions are affected differently 
by the QCD corrections. In this paper, we provide the next-to-leading order (NLO) QCD 
corrections to Higgs pair production including dimension-6 operators in the limit of large 
top quark masses. Depending on the dimension-6 coefficients entering the Lagrangian, the 
new operators affect the relative NLO QCD corrections by several per cent, while modifying 
the cross section by up to an order of magnitude. 
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1 Introduction 


With the discovery of the Higgs boson |1,2], its role has developed from the long-sought particle 
into a tool for exploring beyond the Standard Model (SM) physics |3j, possibly paving the way 
into new physics (NP) territory. This is the more true, as to date we are lacking any direct 
evidence of physics beyond the SM (BSM). The Higgs boson itself behaves SM-like with its 
couplings to other SM particles being close to the predicted values, in particular the couplings 
to gauge bosons. In some NP models, however, the trilinear Higgs self-coupling can still deviate 
significantly from the SM expectations [4|. A means to describe BSM physics, that is realized 
at a scale well above the electroweak symmetry breaking scale, in a rather model-independent 
way is given by the Effective Field Theory (EFT) framework. Deviations from the SM are 
parametrized by higher-dimensional operators, which lead to modifications of the Higgs boson 
couplings to the other SM particles and to itself. 


The trilinear Higgs self-coupling is accessible in double Higgs production 
dominant production mechanism at the LHC provided by gluon fusion |9[10 


with the 
The leading 

order process is mediated by top and bottom quark triangle and box diagrams. As in single 
Higgs production 111,12j, the next-to-leading order (NLO) QCD corrections to this process are 
important. They have first been obtained in the limit of large top quark masses |5j. While 
this approximation works quite well for single Higgs production, the uncertainties of the ap¬ 
proximation are considerably larger for double Higgs production and even more in the case of 



order within the SCET approach. 


Higher-dimensional operators relevant for Higgs pair production through gluon fusion have 
been discussed in [4p3 - 26 . They lead to contributions that are different for the triangle and box 
diagrams mediating the pair production process. In this work we perform the computation of the 
NLO QCD corrections to gluon fusion into Higgs pairs including higher dimensional operators 
in the large top mass limit. Our result allows us to investigate the validity of an approximation 
applied in previous works. This approximation relies on the multiplication of the full leading 
order (LO) cross section by an overall iL-factor, given by the ratio of the SM result for the NLO 
QCD cross section divided by the LO cross section, in the large top mass limit. 

In the next section[2]we present the details of our calculation. This is followed by a numerical 
analysis in section i ^ section [4] we summarize and conclude. 


2 Details of the calculation 

Gluon fusion into Higgs pairs is mediated by top and bottom quark loops dominantly [9|. We 
compute the NLO QCD corrections in the heavy top quark limit and we neglect in the following 
in this framework the bottom quark loops, which only contribute with less than 1% [5 

If physics beyond the SM appears at some high-scale, NP effects can be parametrized in a 
rather model-independent way by introducing higher-dimensional operators. In case the Higgs 
boson is embedded in an S'I7(2)-doublet H the leading BSM effects are given by dimension- 
6 operators^] In the Strongly-Interacting-Light Higgs (SILH) basis the operators relevant for 

1 In certain parameter regions dimension-8 operators can become more important than the dimension-6 ones [ 4 ]. 
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Higgs pair production are given by 27 , 
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( 2 . 1 ) 


where v is the vacuum expectation value v ~ 246 GeV, Mh = 125 GeV the Higgs boson mass, 
mw the W boson mass, yt the top Yukawa coupling constant, g s the strong coupling constant 
and the gluon field strength tensor. Note that we neglect CP-violating effects. An estimate 
of the size of the coefficients ch,c u ,cq and c„ and the most important experimental bounds 


can be found in 28 . The first three operators in Eq. (2.1) modify the top Yukawa and the 


trilinear Higgs self-coupling with respect to the corresponding SM values, while the last operator 
parametrizes effective gluon couplings to one and two Higgs bosons not mediated by SM quark 
loops. The second operator furthermore introduces a novel two-Higgs two-fermion coupling 1291. 


In case the SU(2)l x U(1)y is non-linearly realized and the physical Higgs boson h is a 
singlet of the custodial symmetry and not necessarily part of a weak doublet, the contributions 
relevant for our process are summarized by the non-linear Lagrangian 


30 


A/l n0 n-lin 7) Tflf tt ( Q ^ T Ctt ^2 


h YYc 3 H 3M * 


6 V 


)'> 3 +v g “'“' g ^( c 4 +c 


h 2 \ 

99 YY ) ’ ( * 2 - 2 ) 


' 2v 2 J 


with a s = ^/(47t). In contrast to the SILH parametrization, where the coupling deviations from 
the SM are required to be small, the non-linear Lagrangian is valid for arbitrary values of the 
couplings Q. From the SILH Lagrangian in the unitary gauge and after canonical normalization 
the relations between the SILH coefficients and the non-linear coefficients Cj can be derived, 
leading to |4 


_ 1 CH 

Ct — 1 2 ^ u 5 


Ctt = — - (CH + 3 c u 


C3 = 1 - y C H + C(5 , 


'9 ~ C 99 ~ C 9 ( a ^j ) ( 2 -3) 


with «2 = \/2Gp'rn^y/7i and Gf denoting the Fermi constant. In the following we will give 
results for the non-linear parametrization and defer the SILH case to Appendix |A| 

In the low-energy limit of small Higgs four-momentum an effective Lagrangian valid for light 
Higgs bosons can be derived for the Higgs boson interactions. The effective Lagrangian can 
be used to compute the QCD corrections to Higgs pair production in the large top mass limit. 
From single-Higgs production it is known that the JL-factor obtained in this limit approximates 
the exact value to better than 5%, when the full mass dependence is included in the LO cross 
section 11 . The low-energy approach has also been used to derive the QCD corrections to 
Higgs pair production |5|. Here the uncertainty of 10% induced in the K -factor is worse than in 
the single Higgs case 15-18 The effective Lagrangian for the Higgs couplings to gluons and 

Since the investigation of the concerned kinematic regions is challenging we neglect those operators in the following. 

2 Note, however, that the application of the LET in minimal composite Higgs models at LO leads to an even 

worse approximation than in the SM 14][29). Since the NLO corrections are dominated by soft and collinear gluon 
effects, the top mass effects on the A-factor can be expected to be of order 10-20% also for models including 
higher-dimensional operators. 
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quarks reads |5| 
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(2.4) 


The factor (1 + ll/4a s /7r) emerges from the matching of the effective to the full theory at NLO 
QCD. The Feynman rules for the effective couplings of two gluons to one and two Higgs bosons, 
based on the low-energy theorems (31,32 , are given in Fig. [I] 
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Figure 1: Feynman rules for the effective two-gluon couplings to one and two Higgs bosons in the heavy quark 
limit, including NLO QCD corrections. The incoming four-momenta of the gluons are denoted by hi and fe- 


The generic diagrams contributing to gluon fusion into Higgs pairs at LO are depicted in 
Fig- S The LO partonic cross section can genetically be cast into the form 




|C' A Fi+F 2 | 2 + |c t 2 G' D | 


Jt _ 256(2vr) 3 

where hr denotes the renormalization scale. The Mandelstam variables are given by 


(2.5) 


s = Q 2 


t = Mt - 


2 Q 2 (l — {3 cos 8) 


and 


~ j, j 2 Q 2 {l + Pcos0) 

u = M h --- , (2.6) 


in terms of the scattering angle 8 in the partonic center-of-mass (c.m.) system with the invariant 
Higgs pair mass Q and the relative velocity 


/ 4 M. 2 


The integration limits at cos 8 = ±1 are 


t± = Mt - 


2 Q 2 { l=F/ 3 ) 


(2.7) 


( 2 . 8 ) 
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Figure 2: Generic diagrams contributing to Higgs pair production in gluon fusion at LO. 


The form factors F A and Fa in F\ and F 2 defined as 

2 2 
Fi = c t F A + -ca and F 2 = c\F n + c tt F A - -c D , (2.9) 

contain the full quark mass dependence and can be found in |9j. In the heavy quark limit the 
form factors F A , F n and Go approach 

F a —> — , Fa > — — and Go = 0 , (2-10) 

o o 

and F\ and F 2 simplify to 

= l(c t + c A ) , F^= 2 -(-c 2 t +Ctt-cu). ( 2 . 11 ) 

We have furthermore introduced the abbreviations 

M 2 

Ca = ^hhh Q 2 _ M 2 ^ iMhT] ’ CA = Cg an< ^ Cd = ~^ c 99 ’ (2.12) 


with 


hhh 


M 2 Z 


(2.13) 


The terms proportional to q. respectively c| in F\ and F 2 and in front of the form factor Go are 
the usual SM contributions including the modifications due to the rescaling of the top Yukawa 
coupling by c*. The contributions coming with c A and Co originate from the effective two-gluon 
couplings to one and two Higgs bosons, while the term involving ctt is due to the novel two-Higgs 
two-top quark coupling. 

We use the effective couplings to compute the NLO QCD corrections to Higgs pair produc¬ 
tion. They are composed of the virtual and the real corrections. Sample diagrams are shown 
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Figure 3: Sample effective diagrams contributing to the virtual (upper) and the real (lower) corrections to gluon 
fusion into Higgs pairs. 


in Fig. [3j The calculation is performed by applying dimensional regularization in d = 4 — 2e 
dimensions. The ultraviolet divergences are canceled through the renormalization of the strong 
coupling constant in the MS scheme with five active flavours, i.e. with the top quark decoupled 
from the running of a s . The infrared divergences are canceled by summing the virtual and the 
real corrections. The remaining collinear initial state singularities are absorbed into the NLO 
parton densities, which are defined in the MS scheme with five light quark flavours. The finite 
hadronic NLO cross section can then be written as 


cnlo (jPP ->• hh + X) = (j LO + Acr virt + Aa gg + A a gq + Aa qq . 


(2.14) 


We obtain for the individual contributions of Eq. (2.14) 
OLO 
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where s denotes the hadronic c.m. energy and 


ro = 


4 Ml 


( 2 . 20 ) 


The Altarelli-Parisi splitting functions are given by |33j, 

1 \ 1 


P99 ( z ) — ^ ( 1 
P (,) - 4 1±(1 zz£ 

19qK) ~ 3 z 


H-2 + z( 1 — z) 


33 - 2 N f 

+- q—^S(1-z) 


( 2 . 21 ) 


with Nf = 5 in our case. We denote the factorization scale of the parton-parton luminosities 
d£ 13 /dr by hf- While the relative real corrections are not affected by the higher-dimensional 
operators, the virtual corrections are altered compared to the SM case because of the overall 
coupling modifications of the top Yukawa and the trilinear Higgs self-coupling and due to the 
additional contributions from the novel effective vertices. The coefficient C for the virtual 
corrections reads 


C = 


0 33-2AY u% 11 

” + og (y + ^ 


+ Re 


6 Q 2 2 

dt [CaFi -I- F 2 \ [-C^^ca + ^c n + ci(c t + ca) 2 ] 


Jf + di \CaF 1 + F 2 \ 2 + \bfG a \ 


( 2 . 22 ) 


+c 2 (c t + ca) 2 Re 


f + di c 2 ^(Q* - 2M 2 )G a 
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\C A Fi+F 2 \ 2 + \c 2 t G L 


with the squared transverse momentum 

(t - M 2 )(u - M 2 ) 


Pt = 


and the coefficients 


Q 2 


Cl = c 2 = 


— Mu 
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(2.23) 


(2.24) 


The third line in Eq. (2.22) and the terms proportional to ci in the second line originate from 
the third diagram with the two effective Higgs-two-gluon couplings in Fig. [ 3 ] (upper), and the 
remaining terms are due to the diagrams with gluon loops in Fig. [ 3 ] (upper). In the derivation 
of the coefficient C for the virtual corrections we have kept the full top quark mass dependence 
in the LO amplitude. The SM result is recovered for 


ct 


c t t 


0 , c 3 


ca —1 0 and ca ~> 0 . 


(2.25) 


3 Numerical Analysis 


For the numerical analysis we have implemented the LO and NLO Higgs pair production cross 
sections including the contributions of dimension-6 operators, as presented in the previous sec¬ 
tion, in the Fortran program HPAIR 34 . We have chosen the c.m. energy yfs = 14 TeV and 
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for comparison also the very high energy option y/s = 100 TeV. The Higgs boson mass has 
been set equal to My, = 125 GeV 135] and for the top and bottom quark masses we have cho¬ 
sen mt = 173.2 GeV and mb = 4.75 GeV, respectively. We have adopted the MSTW08 [36] 
parton densities for the LC0 and NLO cross sections with a s (Mz) = 0.13939 at LO and 
a s (M z ) = 0.12018 at NLO. 

In order to study the impact of the new couplings on the QCD corrections we show in Fig. [4] 
the ii-factor, defined as the ratio of the NLO and LO cross sections, A = onlo/olO) where the 
parton densities and the strong couplings a s are taken at NLO and LO, respectively. Deviations 
with respect to the SM A-factor arise in the virtual corrections due to the second and third line 
in the formula Eq. (2.22) for the coefficient C. Additionally, the real corrections are affected 
because of the different weights in the r integration due to the modified LO cross section. In 
Fig. [4] (upper) we have set all couplings but c g to their SM values, C3 = ct = 1 and ctt = c gg = 0. 
We have varied c g away from its SM value c g = 0 in the range —0.15 < c g < 0.15. For illustrative 
purposes we have chosen a rather large range, that goes beyond current experimental limits 
obtained under certain assumptions |4 37 . In the lower plot we have instead set c g = 0 and 
varied c gg away from its SM value c gg = 0 in the range —0.15 < c gg < 0.15 (4j. As can be 
inferred from Fig. [4] the new contact interactions c g and c gg each vary the K -factor between 1.8 
and 1.9 in the investigated range. Figure [4] (lower) shows that the maximal deviation from the 
SM K -factor is obtained for c gg = —0.15, where the A'-factor practically becomes constant. It 
amounts to 


X C 99 
u max 


max] A c »s — A’ sm | 
J^SM 


= 5.4% 


(3.26) 


The impact on the total cross section, however, is much larger. Here we have at NLO max|<r c " — 
(j SM |/cr SM = 5.8. While the effect of higher-dimensional operators on the A-factor is at the level 
of a few per cent, on the cross section itself it is enormous. In Fig. [4] (upper) the maximum 
is less pronounced, as c g modifies the coupling of a single Higgs boson to two gluons, which is 
attached to the Higgs propagator. At the c.m. energy of y/s = 100 TeV (not shown here) the 
A'-factor is smaller and varies for non-zero c gg in the range 1.46-1.58. For c g the impact at very 
high energy is smaller, i.e. a per cent effect. 

In Fig. [5] we have set all couplings to their SM values but allowed for the new contact 
interaction between two Higgs bosons and two top quarks parametrized by ctt, which we have 
varied between -1.5 and 1.5 [4j. The maximum A"-factor is reached for ctt ~ 0.7 where the LO 
cross section is minimized and we have 


= 3.2% . (3.27) 

The maximum deviation is smaller than for c gg . Note, furthermore, that the ctt value, for which 
the deviation is maximal, is much larger than for the c g and c gg variations discussed before. 
This is due to the normalization of these coupling factors as can be inferred from the Feynman 
rules in Fig. [lj At y/s = 100 TeV the A-factor varies between 1.49 and 1.59. 

The value of the trilinear coupling is practically not constrained and we have allowed in 
Fig. § for a variation of C 3 between -10 and 10 |4|, while setting all other couplings to their SM 

3 Note that in the LO cross section we have neglected the bottom quark loops. Their effect amounts to less 
than 1%. This treatment accounts for the fact that we did not introduce a coupling modification factor in the 
Higgs Yukawa coupling to bottom quarks. Additionally it is consistent with the application of the heavy top 
quark limit in the NLO QCD corrections. 
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Figure 4 : IL-factors of the QCD-corrected gluon fusion cross section cr(pp —» hh + X) at the LHC with c.m. energy 
^/s = 14 TeV. The dashed lines show the individual contributions of the four terms of the QCD corrections given 
in Eq. ( 2 . 141 , i.e. Ki = Acri/oLO (* = virt,gg,gq,qq). We have set the SM values C 3 = Ct = 1, Ctt = 0 and varied 
c g with c gg = 0 (upper), respectively, varied c gg with c g = 0 (lower). 
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values. The impact of C 3 is rather small. In the investigated range the maximum deviation of 


5 C3 

max 


= 2 . 1 %, 


(3.28) 


is reached for 5 < C 3 < 10. At ->/s = 100 TeV the AT-factor varies between about 1.42 and 1.51. 
Finally, the variation of q in the still allowed range 0.65 < a < 1.15 37 , while keeping all 


other couplings at their SM values, hardly changes the A'-factor, and we do not show the corre¬ 
sponding results separately. In this parameter configuration ct enters in the virtual correction 
factor C, Eq. ( 2 . 22 ), through the terms proportional to ci and C 2 , and its effect in the nominator 
almost cancels against the one in the denominator. 


The above discussion shows, that the new couplings affect the A'-factor by only a few per 
cent. It has to be kept in mind though that we varied the couplings only one by one away from 
their SM values. The combination of various new couplings could have a larger impact on the 
NLO corrections. 


4 Conclusions 

We have computed the NLO QCD corrections to Higgs pair production including dimension- 
6 operators in the large top mass limit. The dimension -6 operators lead not only to coupling 
modifications of the SM Higgs couplings, but also induce effective gluon couplings to one and two 
Higgs bosons and a novel two-Higgs two-top quark coupling. The various contributions to Higgs 
pair production are affected differently by the QCD corrections. Depending on the relative 
size of the NP contributions, the A'-factor is changed by several per cent in the parameter 
regions compatible with the LHC Higgs data. This small impact on the A"-factor underlines the 
dominance of soft and collinear gluon effects in the QCD corrections. The inclusion of the QCD 
corrections in the gluon fusion process based on the effective theory approach to describe NP, 
is necessary for reliable predictions of the cross section. 


Appendix 

A Gluon Fusion into Higgs Pairs in the SILH Approximation 

The SILH approximation of NP effects is valid for small shifts bc% in the Higgs couplings c* away 
from the SM values c? M , i.e. 

SILH: Ci = cf M (l + 6Ci) , with §Ci <C 1 . (A.29) 

While in the non-linear case arbitrary values for the coupling coefficients are allowed and terms 
quadratic in bc{ have to be taken into account in order to avoid non-physical observables, such 
as negative cross sections, in the SILH approach an expansion linear in 6ci has to be performed. 
With 

CH + 2c u _ CH + 3 c. u 

ct — IT oct — 1 ^ > ctt — octt — 2 j 

c 3 = 1 T 8c 3 = 1 — 3Cg 9 2CS , Cg = 5c g = 5c gg = Cg 3 (A.30) 
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cf. Eq. (2.3), this yields for the LO partonic cross section Eq. (2.5) in the SILH parametrization 
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The NLO SILH cross section is obtained from Eqs. (2.14)-(2.19) by replacing 
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